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A new model for chemical reactions, the Quantum-Kinetic (Q-K) model of 
Bird, has recently been introduced that does not depend on macroscopic rate 
equations or values of local flow field data. Subsequently, the Q-K model has 
been extended to include reactions involving charged species and electronic 
energy level transitions. Although this is a phenomenological model, it has 
been shown to accurately reproduce both equilibrium and non-equilibrium 
reaction rates. The usefulness of this model becomes clear as local flow 
conditions either exceed the conditions used to build previous models or 
when they depart from an equilibrium distribution. Presently, the 
applicability of the relaxation technique is investigated for the vibrational 
internal energy mode. The Forced Harmonic Oscillator (FHO) theory for 
vibrational energy level transitions is combined with the Q-K energy level 
transition model to accurately reproduce energy level transitions at a 
reduced computational cost compared to the older FHO models. 


Nomenclature 

c r = Relative collision velocity (m/s) 

E = Energy of particle (J) 

/ = Distribution function 

g = Energy level degeneracy 

g* = Constant in vibrational transition probability (m/s) 

k = Boltzmann’s constant (1.380 650 x 10' 23 J/K) 

m r = Reduced mass in collision (kg) 

N = Number of particles 

P = Total probability of transition 

T = Temperature (K) 

8 = Dirac delta function 

e = Average number of quanta transferred between vibrational levels 

r(x) = Gamma function of x 

n , J = Probability of transitioning from level i to j 
0 = Characteristic temperature (K) 

co = Viscosity index for VHS potential 
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Subscripts 


coll 

= Collision energy 

int 

= Internal energy 

i,j,k,m 

= Internal energy level level indices 

max 

= Maximum internal energy level possible 

tr 

= Translational energy 

trans 

= Transition 

vib 

= Vibrational energy 


I. Introduction 

Tntemal energy relaxation is an important physical process in hypersonic flows. While the 
-^-rotational energy mode is a concern in all flow regimes, the vibrational and electronic energy 
modes become more significantly activated in the high-temperature flows generally associated 
with hypersonic flows. To complicate matters, the relaxation rate of the internal modes is finite, 
which may lead to a state of thermal non-equilibrium between the translational, rotational, 
vibrational, and electronic temperatures. Under these conditions, it is also possible for the 
distribution of the quantized internal energy modes to not follow the equilibrium Boltzmann 
form. 

The method of computation of hypersonic flows depends on where along the trajectory of the 
vehicle the solution is desired. When the vehicle is in the continuum flow regime, computational 
fluid dynamics (CFD) is used. For non-continuum flow (i.e. rarefied to transitional), the direct 
simulation Monte Carlo (DSMC) method is generally used because of the breakdown of the 
continuum assumptions and thus kinetic methods become more appropriate. In contrast to CFD, 
DSMC performs a direct physical simulation of the gas at the molecular level. In the simulation, 
particles are tracked in space and time, accounting for both gas-surface interactions and 
intermolecular collisions. 

Finite rate rotational and vibrational energy mode relaxation is usually modeled in the DSMC 
method using a probability of energy exchange for each collision that is consistent at the 
macroscopic level with the derivations of Parker[l] and Landau and Teller[2], respectively. For 
the fraction of collisions that are selected for rotational or vibrational relaxation, a quantized 
form[3, 4] of the Larsen-Borgnakke[5] (LB) energy exchange model is employed. These 
approaches act to equilibrate the translational and rotational/vibrational energy modes of the 
molecules without regard to transition rates between the individual quantum levels. 

Beyond the simple models mentioned above, there has also been work towards state-to- state 
energy exchange models for both rotational [6-9] and vibrational [10- 12] relaxation. These 
models make use of more detailed evaluations, such as quasi-classical trajectory analysis or 
forced harmonic oscillator models, to calculate transition rates between quantum levels as a 
function of the relative velocity/energy of the colliding particles in an attempt to more accurately 
reproduce the non-equilibrium distribution of the energy levels. These high-fidelity models have 
allowed much better comparison to measured transition rates and energy level distributions. 

Recent advances in DSMC algorithms [13- 15], specifically the Quantum-Kinetic (Q-K) 
chemistry model, have allowed the prediction of equilibrium and non-equilibrium reaction rates 
using only kinetic theory and fundamental molecular properties (i.e., no macroscopic reaction- 
rate information is needed). Models such as this are important especially when reaction rate data 
are limited or the gas is far from equilibrium and an engineering approximation is required. 
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Following the Q-K model, an electronic energy relaxation model was introduced[16-18] that was 
shown to reproduce energy level transition rates within the range of transition rates reported in 
the literature. This model was subsequently added to NASA’s DSMC code DAC (DSMC 
Analysis Code) [19] and compared to CFD predictions of the Fire II flight test[20]. It is the 
purpose of this study to investigate the applicability of the relaxation model developed for 
electronic energy to the vibrational energy mode. 

II. Extending the Quantum-Kinetic Chemistry Model 

One of the attractive properties of the Q-K model of Bird is that it does not rely on any 
macroscopic information, leveraging instead off of kinetic theory and physical data 
corresponding to the particles undergoing each individual collision. The new methods that were 
developed for treating electronic energy level transitions will be generalized to any internal 
energy mode being described by quantum levels as originally derived. 


A. Detailed internal energy level model 

In order to implement a model for the distribution of internal energy in the DSMC technique, 
there are three procedures that must be defined. First, when a particle is introduced into a 
simulation, it is necessary to obtain a new energy level from a given distribution through 
statistical sampling. A new energy level is also required when a particle changes its energy after 
a collision with a surface. Under equilibrium conditions, the distribution has the well-known 
Boltzmann form. Second, it is necessary to statistically sample a new energy level following a 
collision, which involves internal energy transfer. Third, a method to reproduce the internal 
energy transition rate must be defined. 

1. Equilibrium sampling 

Each energy level j has a distinct energy, Ej, and degeneracy, gj. The Boltzmann distribution 
for the internal energy levels at a given temperature T gives the following result for the fraction 
of particles in level j : 


= N j= 8j exp(-Ej/kT) 

^' nt AT 'max 

£& ex v(-EjkT) 


( 1 ) 


This distribution is used when either creating a new particle or after a surface collision in the 
DSMC simulation at a boundary specified at the temperature T. However, it is not possible to 
sample an internal energy level j directly from the distribution. Therefore, an acceptance- 
rejection procedure is used. This is performed by selecting values for j from the following 
distribution: 


J int \J) ( ~W I / 7 m \ 

g, exp(- E, / kT ) 


( 2 ) 


where J is the value of j with Equation (1) maximized. The sampling of a new internal energy 
level then proceeds, in general, as follows: 

1 . select at random an energy level, j, evenly distributed between 0 and J max , where J max 
is the maximum possible energy level; 

2. determine the value of J with the Boltzmann distribution in Equation (1) a maximum; 
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3. accept the value of j if f' nt ( j) > RAND , a random number 

4. if the value of j is not accepted, return to step 1. 

2. Post-collision sampling 

A phenomenological approach is usually adopted in the DSMC method when a collision 
occurs that involves energy transfer. The LB method samples a post-collision state from a 
combined distribution of the translational and internal collision energies of the colliding 
particles. Based on the approach of Bemgemann and Boyd[4] and Boyd[3], the Dirac delta 
function is used to determine the distribution of energies in Equation (1) in the following form: 




£& txp(-E,/kT) 


{ kT 


int 

kT 


;=o 


(3) 


Consideration must now be given to the distribution of translational energy of the colliding 
particles. This distribution is naturally affected by the intermolecular model used since this 
determines the collision probability. For the present study, the variable hard sphere (VHS) 
collision model of Bird[21] is used. However, it is a simple matter to develop the formulation 
for an alternative collision model such as the variable soft sphere model of Koura[22], For the 
VHS model, the distribution of the translational energies is: 


/„(£„)= 


r (%-o>) 


f E \ 
kT J 


exp 


E \ 

n tr 


kT 


(4) 


Using Equations (3) and (4), the combined distribution for sampling a post-collision internal 
energy level j’ from the total collision energy E coll = E tr + E mt = E' r +E' nt is: 


f(f’ E cou) = 


Si 
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^ coll 


■ F J 

^int 


X Si exp (- E i /kT ) L 


kT 


exp 
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^ coll 
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(5) 


In applying the general LB scheme, it is assumed that local thermodynamic equilibrium 
prevails. Therefore, the temperature T in Equation (5) is constant. Also, the total collision 
energy is constant, so it is only necessary to perform sampling of the post-collision state for j’ 
from the following distribution form: 

s(j'-.E,A^g r (E„,-ELf~" ( 6 ) 


Again, an acceptance-rejection procedure is used. The normalized distribution that is required is 
obtained by finding the value of j’ for which Equation (6) is a maximum, J’, which can be 
different for each value of E co u. Therefore, the following distribution is obtained: 




S/ ( E co u ^Lt) 
Sr (Ecou ~ E ml ) 


( 7 ) 


In the DSMC code, these procedures are implemented in the following way: 
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1. given a pair of particles with total collision energy E co u, which undergoes internal 
energy exchange, determine J’ from Equation (6); 

2. determine the maximum allowable internal energy level obtainable from E co u, J max \ 

3. take J* to be the smaller of J’ and J max ', 

4. as described in the previous section on equilibrium sampling, perform an acceptance- 
rejection procedure to sample j’ from 




S f [Ecou 


( 8 ) 


3. Internal energy level transitions 

A transition from internal energy level i to level j for species A can be written in the form of a 
chemical reation as A 1 + Z? — > A 1 + B . Since we are treating the transition as a chemical reaction, 
there are J max transitions to consider where J max is the maximum energy level possible. The total 
probability of the transition occurring is then: 


pij 

1 trans 


= n 


Sj(e co 11 -E:; m ) 


int J m 


I g,( E „.- E LY 


k=0 


(9) 


where n‘; J t represents the probability of transitioning from level i to level j (discussed 
subsequently), and the last term is the probability of choosing level j. In the DSMC code, this 
process is broken into two parts. First, the energy level transition is chosen by the method 
outlined in Section 2 above for post-collision sampling. Then, the possibility of transition to that 
state is tested by comparing the probability fl';^ to a random number. The process for 
determining if a transition is to occur and what level to transition to is then as follows: 

1. determine the value of J max and compute the denominator of Equation (8); 

2. perform acceptance-rejection according to Equation (8); 

3. calculate the probability II ' n ( ; 

4. if the probability is greater than a random number between zero and one, the transition 
occurs. 


B. Application of internal energy model to specific modes 

Now that the general procedure for internal energy mode relaxation has been outlined, the 
procedure is applied to the electronic and vibrational internal energy modes. The new internal 
energy level transition model as discussed above is investigated by its ability to reproduce 
transition rates from the literature. The DSMC models are applied in a three-dimensional code 
where the results are accumulated over all of the cells (1,000 particles per cell in 1,000 cells) to 
approximate a zero-dimensional simulation. The test gas is molecular nitrogen and temperatures 
ranged between 2,000 K and 20,000 K at a number density of 1.0e21/m 3 in an adiabatic cube 
measuring 1mm on a side. It should be noted that the uncertainty of the currenly available 
measured/calculated rates can be very large. The uncertainty is not always reported and in some 
cases is unknown. Generally, if more than one source is available for any given rate, the spread 
is as much as an order of magnitude or, for many cases, more. 
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1. Electronic energy levels 

Very little simplification can be made for the 
electronic energy level transition model. The 
electronic energy mode has unequally spaced 
levels that each have their own degeneracy. This 
energy mode was the topic of earlier studies [lb- 
18] and the model described herein is nearly as it 
was in the references. In the electronic energy 
mode, the probability of transitioning, Yl l f le , was 
originally set by selecting a random possible 
energy level and then if the randomly selected 
energy level was equal to the level selected by 
step (2) above, the transition was allowed to 
occur. This is equivalent to setting the transition 
probability equal to the inverse of the number of 
possible energy levels. The resulting sample 
transition rate was found to compare well with 
rates from the literature. An example is 
presented in Figure 1 for an N 2 -N 2 collision for a 
variety of energy level transitions (0-1, 1-2, 2-3) 
where the sampled DSMC transition rates are compared to data from Park[23], Flagan[24] and 
Nagy [25]. Note that the spread in reported transition rates can be quite large. More examples, 
including other collision partners, are presented in References [16-18]. 



Temperature (K) 

Comparison of sampled and 
reported electronic energy level transition 
rates for collisions of molecular nitrogen. 


Figure 1. 


2. Vibrational energy levels 

The vibrational energy levels are usually assumed to be described by a simple harmonic 
oscillator (SHO). While this does not have to be the case, it is overwhelmingly chosen because it 
is simple to model. However, the SHO model for vibrational energy assumes that the energy 
levels have a degeneracy of one and are equally spaced by the energy defined by A E vib = Q vib k , 

such that the energy in level j is defined as E J vib = j® vib k . Therefore, the total probability of 
transition from Equation (9) becomes: 


pij _ rr' j _ 

1 vib 1 l vib J, 


i^coll jk®vib) 


'Z( E c«ll- mke ,b) 


y 2 -a 


m = 0 


( 10 ) 


Comparisons of the sampled DSMC vibrational energy level transition rates are made to rates 
quoted by Hansen[26] for excitation and by Billing[27] for de-excitation. When the probablility 
is set to the inverse of the number of available states as was done for electronic energy, the 
resulting transition rates are much too fast as presented in Figures 2 and 3 (Q-K: Original). 
However, when the Millikan- White[28] formulation of Bird[29] is used, the resulting transition 
rates are slightly too slow, as has been found by various researchers [10, 12], 
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Figure 2: Comparison of sampled and reported vibrational energy level transition rates for 
collisions of molecular nitrogen between levels 0 and 1 (left) and 1 and 0 (right). 
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Figure 3: Comparison of sampled and reported vibrational energy level transition rates for 
collisions of molecular nitrogen between levels 2 and 3 (left) and 3 and 2 (right). 


The probability of transitioning from level i to level j for the vibrational mode can be found 
from sources such as the small perturbation theory of Hansen [26] or the forced harmonic 
oscillator (FHO) model of Adamovich[30], presented for the DSMC method by 
Vijayakumar[31]. For example, the probability of transitioning from level i to level j according 
to the FHO model is: 


n'i = /!./!e^exp(-e) 


f (~l) r 1 2 

Ao r\(i-r)\(j-r) \ £ r 


( 11 ) 


and 


7 

American Institute of Aeronautics and Astronautics 


( 12 ) 


£ = 
a = 




-g*!c r 


m r g 


*2 


ke 


vib 


In these two equations, m r is the reduced mass of the colliding pair, k is Boltzmann’s constant, 
Q vib is the characteristic vibrational temperature of the relaxing molecule, c r is the relative 
velocity of the colliding pair, and g* is a constant suggested as 62,000 m/s in Ref. [31]. When 
this probability of transition is used, better comparison is found to the rates from the literature as 
presented in Figures 2 and 3 (Q-K: FHO). 

The advantage of using the Q-K model presented herein over the original model presented by 
Vijayakumar is computational speed. In the original model, the transition probabilities of each 
energy level energetically possible must be computed, which could be up to 34 levels for 
nitrogen. Instead, the current model requires only one of the probabilities to be computed, so the 
speedup is of the order of 10 times the original. The computational time is comparable to the 
Millikan- White relaxation model, but comparisons to data from the literature are much better. 


III. Conclusion 

Extensions to the Quantum-Kinetic chemistry method for direct simulation Monte Carlo 
simulations have been proposed for internal energy state-to-state transitions. Comparisons have 
been made for electronic energy level transitions where the proposed transition rates are within 
the scatter of the reported rates in the literature and for vibrational energy level transitions where 
the proposed transition rates compare well with reported rates. The primary advantage of the 
current transition rate prediction model over previous models is, as was the original Q-K model, 
the simplicity of implementation and the need for only a few basic properties of the particle with 
the attendant computational speedup. 

Several improvements to the proposed model are needed. First, rotational energy should be 
included in the model. The equilibrium and post-collision energy state assignment routine 
outlined above is applicable to the quantum model for rotation, but an appropriate energy level 
transition probability model has not yet been identified. As with the vibrational energy model, 
an approximation, such as Boyd’s collision energy based transition model[3], could be used in 
place of a state-to-state transition model. The second improvement to the proposed model would 
be a more rigorous treatment of the electronic energy level transition model by making use of 
state-to-state transition probabilities such as that of Carlson[32], However, using a methodology 
such as this is counter to the simplicity of the Q-K model in that specific reaction rate data would 
be needed. 
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